Multibranch entrainment and slow evolution among branches in 

coupled oscillators 

Toru Aonishi 1 and Masato Okada 2 ' 1 

1 Brain Science Institute, RIKEN, 2-1 Hirosawa, 
Wako-shi, Saitama, 351-0198, Japan 
2 ERATO Kawato Dynamic Brain Project, 2-2 Hikaridai, 
Seika- cho, Soraku-gun, Kyoto 619-0288, Japan 
(Dated: February 1, 2008) 

Abstract 

In globally coupled oscillators, it is believed that strong higher harmonics of coupling functions 
are essential for multibranch entrainment (MBE), in which there exist many stable states, whose 
number scales as ~ O(expiV) (where N is the system size). The existence of MBE implies the 
non-ergodicity of the system. Then, because this apparent breaking of ergodicity is caused by 
microscopic energy barriers, this seems to be in conflict with a basic principle of statistical physics. 
In this paper, using macroscopic dynamical theories, we demonstrate that there is no such ergodicity 
breaking, and such a system slowly evolves among branch states, jumping over microscopic energy 
barriers due to the influence of thermal noise. This phenomenon can be regarded as an example of 
slow dynamics driven by a perturbation along a neutrally stable manifold consisting of an infinite 
number of branch states. 

PACS numbers: 05.45.Xt, 75.10.Nr, 87.18.Sn 
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An Ising model with a ferromagnetic interaction has a simple phase space, and below 
some critical temperature it only possesses two macroscopic states, two oppositely directed 
ferromagnetic states which have the same energy. Contrastingly, in the case of ferromagnetic 
coupled oscillators with coupling functions containing strong higher harmonics, each oscilla- 
tor is multi-stable, and as a result, the configuration observed in a simulation will consist of 
individual clusters characterized by different phases. Daido showed through numerical sim- 
ulations that there exist many stable states, whose number scales as ~ O(expiV) (where N 
is the system size)0. This complex mutual entrainment is called multibranch entrainment 
(MBE). This phenomenon share the same mechanism of so-called clustering, which has been 
found in various systems, e.g., globally coupled maps || and pulse-coupled oscillators 0. 
When such systems exhibit periodic behavior, their dynamics can be described by phase 
equations with coupling functions containing strong higher harmonics . 

It is a general principle of statistical mechanics that ergodicity can be broken only by 
macroscopic barriers of the free energy whose heights scale as ~ 0(N). Here, we evaluate 
the free energy of a system of phase oscillators with coupling functions containing strong 
higher harmonics. Then, we demonstrate that MBE means ergodicity breaking caused by 
microscopic free energy barriers whose heights scale as ~ 0(1). As this result reveals, MBE 
seems to be in conflict with a basic principle of statistical physics, because a system can jump 
over these microscopic barriers due to the influence of thermal noise. Next, we apply the 
moment approach [|], [5] to such a system and study properties of its relaxation process. Then 
we elucidate the effect of microscopic energy barriers and demonstrate thermodynamical 
instability of MBE. 

The ferromagnetic Ising model and the ferromagnetic analog network, governed by a 
Markovian process and a continuous process, are in the general class of time- dependent 
Ginzburg-Landau (TDGL) models 0, because in such systems, the dynamics of macroscopic 
order parameters are described by the steepest-descent of the free energy. A ferromagnetic 
coupled oscillator, however, is not in the class of TDGL models. In this paper, we construct 
a Markovian system that has the same free energy as a ferromagnetic coupled oscillator but, 
nonetheless, belongs to the class of TDGL models. Then, comparing this TDGL system and 
non-TDGL ferromagnetic coupled oscillator system, we make clear that microscopic energy 
barriers have no effect on the relaxation of the system in TDGL models we consider but that 
they cause the relaxation to become slow in non-TDGL models we consider. Our results 
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suggest that such slowing of relaxation is characteristic only of non-TDGL models. 

We propose the following phase model with uniform all-to-all coupling by modifying 
Daido's model 



Here fa is the phase of the ith oscillator (with a total of N), and U{ represents its quenched 
natural frequency. The natural frequencies are randomly distributed with a density repre- 
sented by g(co) = (27rcr 2 ) _1 / 2 exp(— co» 2 /2o" 2 ). The parameter a denotes the degree of mixture 
of the fundamental harmonic component and the higher harmonic component in the cou- 
pling functions (0 < a < 1). This model belongs to the class of ferromagnetic systems and 
does not exhibit frustration in the mutual interaction. The function represents a Gaussian 
independent white-noise process satisfying (n,i(t)) = 0, {rii(t)nj(t')) = 2T5ij5(t — t'), where 
T is the temperature of the system, and the inverse temperature is defined as (3 = 1/T. 
Daido used a coupling of the form sin x + a cos 2x instead of that in Eq. (JJ) . In the case of 
Daido's model, when a gradually increase, another stable solution suddenly appears apart 
form original stable solution in coupling function. This sudden appearance of another sta- 
ble solution also occurs in our model. Thus, the bifurcation structure of Daido's model is 
qualitatively equal to that of our model. Moreover, when there is no frequency disorder, 
detailed balance holds in our model. Therefore, we can analyze the MBE of our model in a 
transparent manner. 

For convenience, we define s, = exp(i<fti) and rewrite Eq. ([I]) in the following form with 
bare mean fields: 



Here the overline denotes the complex conjugate. The quantity mi is a mean field consisting 
of the fundamental harmonics {sj}j=i,'-,N, and rri2 is a mean field consisting of the higher 
harmonics {s]} j= i... :N : mi = s j'> m i = jfY.f s r 




3=1 

+a sin (2 (fa - fa) )] + n;(£) . 



(1) 



dt 



Vi + —[(1 - a) (siirii - s,mi) 
+a (s 2 m 2 - sfrfhi)] + n^t). 



(2) 
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Equation (0) obeys the Fokker-Planck equation 

dp(<j)\u},t) d 2 p((p\uj,t) 9 

v(4>, uj) = uj 

+ ((1 - a) (sm 1 - sr%) + a (s 2 m 2 - s 2 m^)) , (4) 

where p(4>\uj,t) denotes the one-oscillator probability density at time t. The site index i of 
uj and can be omitted, because Eq. @j has already been reduced to a one-body system. 

We derive the macroscopic dynamics of the system by using the moment approach j|, [5| . 
We define the moments of p(<ft\uj,t) as 

/•2tt 

X n (uj,t) = / d(f)p((f)\uj , t) exp(— iri(f)). (5) 
Jo 

The motion of these moments obeys the simple equation 

d\ n (uj,t) 



dt 
n(l 



- —Tn X n (uj,t) + inuj\ n (uj,t) 

(mi\ n+ i(uj,t) -mi\ n -x{(jo,t)) 
— —{m 2 Xn+2{uJ,t) -m^X n _ 2 {uj,t)) , (6) 

which is derived through a moment expansion of Eq. (|3]). Equation (^j) possesses a hierarchal 
structure, because the nth moment depends on neighboring moments [the (n— 2)th, (n— l)th, 
nth, (n+ l)th and (n + 2)th moments]. The Oth moment is found trivially to be Ao(^, t) = 1. 
The mean fields ni\ and m 2 are given by rn{ = J dug(uj)Xi (uj, t) , = J dujg(uj)X 2 (uj, t). 
When there is no frequency disorder, we obtain n¥j~ = Ai (0, t) and = X 2 (0, t). If T ^ 0, 
higher moments rapidly decay, since the coefficient of the damping factor in Eq. is 
Tn 2 . Therefore, even if higher moments are omitted in the numerical calculation of Eq. (^|), 
highly accurate results can be obtained. For the convenience of later discussions, we call 
the system without the thermal noise term the "Liouville part", because Eq. (^) without 
the diffusion part (i.e., the damping factor) is equivalent to the so-called Liouville equation. 
The diffusion part in Eq. @ can be regarded as a perturbation of the Liouville equation. 

When there is no frequency disorder [i.e. when g(u) = 6 (no)], detailed balance holds, 
since the coupling function is an odd function |5], 0. Thus, in the limit a — > 0, this model 
can be mapped to the equilibrium system. In this case, there exists a free energy, which 
enables us to apply equilibrium statistical mechanical methods to the system. Then we can 



obtain the Gibbs distribution and partition function 



P(0) = e-Wft/Z, Z = r dcfte-^, (7) 

Jo 

in equilibrium [EL where H denotes the Hamiltonian expressed by 

H(4>) = 

i ^ ' 

We thus obtain the free energy / = — l/(/3N) log Z. Using the saddle-point method to 
evaluate the integral in Z, for large N, we derive the free energy as 

r ( \ 1 — a . 2 . a I 1 2 

f(m 1 ,m2,m 1 ,m 2 ) = — — \mi\ + -\m 2 \ 

_ll g/ e / 3 (^(^l+sW)+f(s 2 m 2 +s 2 m2))\ _ /g-j 

In the limit T — > 0, this becomes 

/(mi, 77i2, ml, ^2) = ^Y~' mi ' 2 + J m2 l 2 = ( 10 ) 

This is the Hamiltonian normalized with respect to N. 

Let us consider the fixed points of Eq. (||) in the limit T — > and a — >• 0. Figure |l| 
displays a typical example of a fundamental harmonic component and a higher harmonic 
component of Eq. (0) in the case ^ = 0. Namely, figure [I] shows two curves hx(<pi) = 
(^ m i — s i m T) an d h 2 ((t>i) = (Si m 2 ~ sfrn^). If a is sufficiently large, Eq. (Q) admits 
four solutions corresponding to intersection points of the fundamental harmonic component 
and the higher harmonic component, as shown in Fig. [I]. The filled circles and open 
circles correspond to stable solutions and unstable solutions, respectively. We thus find that 
this system is microscopically bistable. The two enclosed areas Ai and A 2 in Fig. [I] are, 
respectively, equal to the energies of the two stable solutions, as measured with respect to the 
unstable solution between them. The microscopic energy barrier discussed in introduction is 
the energy barrier between these stable solutions. It is a principle of statistical physics that 
in a situation like this, only the stable solution with the larger enclosed area can be realized, 
because only the stable solution with the lowest energy is thermo dynamically stable. This 
rule is referred to as Maxwell's rule ||. Note that two stable solutions have the same energy 
only if a — 1. However, if system is non-ergodic, multiple states can coexist, and Maxwell's 



rule does not hold. In this case, depending on the initial conditions, each oscillator comes 
to belong to one of two clusters. Then the total number of stable solutions in the whole 
system is 0(2 ) ~ O(expiV). This is the situation involving MBE. 
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FIG. 1: Sketch of microscopic bistability. 




FIG. 2: Macroscopic landscape of the free energy. 

Here, we demonstrate that microscopic multi-valley of the free energy is not reflected in 
the macroscopic structure of the free energy. Figure displays the macroscopic landscape 
of the free energy (Eq. (|9])) as a function of \mi\ and \m2\- The left figure corresponds to 
the case of a monostable system (a = 0.2), and the right figure corresponds to the case of a 
bistable system (a = 0.8). There are no any local minima of the two energy surfaces. Thus, 
it is seen even in the case of a bistable system, the free energy does not possess a multi-valley 
structure macroscopically. From Eq. fli~0l), it is clear that, even in the limit T — > 0, there 
is no macroscopic multi-valley structure in the free energy. Note that in the limit T — > 0, 
the critical value of a, representing the boundary between monostability and bistability, is 
given by a c = 0.467. 

As the macroscopic structure of the free energy reveals, it seems that the ergodicity 
breaking by microscopic barriers does not occur. However, when g(uj) = T — > and 

a > a c , A„(0, oo) = r + (1 — r)(— l) n corresponding to some branch state is a stable solution 
of Eq. (|6p. In this case, the mean fields mi and m2 are given by mi = 2r — 1 and m 2 = 1. 
These moments yield the probability density function p(<f)\0, oo) = r5(<p) + (1 — r)5(<p — 7r), 
where r denotes the ratio of oscillators at <fi = to those at = it. This solution is satisfied 
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for any r (0 < r < 1). In this case, these solutions form a line attractor parameterized by r: 
M = {r + (1 — r)(— l) n |0 < r < 1} . Within this attractor, all of the solutions are neutrally 
stable. We thus find that in the limit N — > oo, the set of solutions forming MBE can be 
parameterized by the single variable r. 

Therefore, in the limit T — > 0, the global ergodicity demonstrated by the macroscopic 
structure of the free energy seems to be in conflict with the existence of M revealed by 
the moment approach. In the following numerical calculations, we evaluate the dynamical 
properties of this system at low temperature, and then answer the question what happens 
in the limit T — ► 0. In the numerical simulations we now discuss, the initial conditions 
of the system were assigned as random numbers with the density function p(<p\uj, 0) = 
qS((p) + (1 — q)5((j) — 7r), where q denotes the ratio of oscillators at = to those at (f> = n. 
Thus, initial conditions were chosen on the manifold M. 
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FIG. 3: Temporal behavior of rrt\ — and ni2 — m^°: (a) a = 0. (b) a = 0.2. The dashed curves 
were obtained analytically, while the solid curves were obtained from numerical simulations. 



Figures |3](a) displays the temporal behavior of \mi — m^°| and \rri2 — iri??\, where mf 3 and 
m|° denote values of m\ and m 2 i n a n equilibrium state corresponding to the minimum point 
of the free energy shown in Fig. |^. The dashed curves in these figures are the analytical 
results obtained from Eq. (|B|), and the solid curves are numerical results. As these figures 
reveal, the analytical results fits in with the simulation results. 

We now calculate the time constant r of the system using logarithm plots. According to 
Fig. H|(a), the relaxation process to equilibrium states is characterized by the exponential 
convergence. The time constant of this convergence, r, for the bistable system (a = 0.8) is 
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about 100 times larger than that for the monostable system (a = 0.2). We also evaluate the 
time constant r as a function of T in the case of a bistable system (a = 0.8). As Fig. |^(a) 
reveals, r diverges to oo in the limit T — > 0. 

Note that in the initial stages of the relaxation process, when initial states are far from 
M, the speed of convergence for the system with a = 0.8 is almost the same as that for the 
system with a = 0.2. During the approach to the equilibrium state, the relaxation process 
suddenly switches from a process of fast dynamics to a process of slow dynamics (see Fig. f|). 
The mechanism responsible for this switch is the following. Since the thermal perturbation 
(the diffusion part in Eq. (|6]) ) is very small, for solutions far from M, the perturbation is 
negligible, and they are attracted rapidly to M, just as in the unperturbed system. (This 
rapid attraction is called the "fast dynamics".) Now note that the Liouville part vanishes 
on M. For this reason, once solutions enter the neighborhood of M, where the effects of the 
Liouville part and the perturbation become comparable, the solution begins to move slowly 
along M toward the ground state. (This is called the "slow dynamics".) 

fast dynamics 

V_ slow dynamics 

r=0 M r=l 



FIG. 4: Sketch of dynamical flow around M. 



0.0(1 
0.000 



(4 



T o.i 



0.001 
0.0001 



(bf 



(c) 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 

T 



FIG. 5: Time constant r as a function of T: (a) a = 0, a = 0.8, q = 0.3; (b) a = 0.1, a = 0.8, q = 0.3; 
(c) Markovian system (a = 0, a = 0.8, q = 0.3). 



In order to compare the non-TDGL system considered above to a TDGL system, we 
evaluated the time constant r for a Glauber system. This system was constructed to have 
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the same Gibbs distribution as Eq. (0) in the equilibrium state. Thus, its macroscopic 
dynamics are described by the steepest-descent of the free energy given in Eq. (|9|). As Fig. 
|5|(c) reveals, in this case r converges to some finite value as T — > 0, even for a bistable system 
(a = 0.8). In the limit T — > 0, the Glauber system chooses only the stable solution with 
the lower energy. We thus find that microscopic energy barriers do not affect the relaxation 
process of the Glauber system. This contrasts with what we have found for the relaxation 
process of the continuous system Eq. ([[]), which is hampered by the microscopic energy 
barriers, resulting in a slow (infinitely slow in the T — > limit) relaxation to the ground 
state. 

In order to study the effect of the frequency disorder, we also measured the temporal 
behavior of \rrii — mf\ and \rri2 — mJf | in the case a 7^ (see Fig. 0(b).) Here, mf and 
denote values of mi and ni2 in a macroscopic steady state [I]. Even if there is frequency 
disorder, the time constant of this convergence, r, for the bistable system (a = 0.8) is about 
100 times larger than that for the monostable system (a = 0.2). As shown in Fig. [|(b), r 
diverges to 00 in the limit T — > 0. 

In conclusion, the system we have considered evolves slowly through a series of branch 
states to the ground state, jumping over microscopic energy barriers through the influence 
of thermal noise. The time constant r characterizing this relaxation process diverges to 00 
in the limit T — > 0. From the macroscopic viewpoint, this phenomenon can be regarded as 
slow dynamics driven by a weak thermal perturbation along the neutrally stabile manifold 
M consisting of an infinite number of branch states. From these results, we can conclude 
that, in the long term, MBE is unstable thermodynamically, but in the short term, the 
system preserves phase patterns highly correlated with the initial conditions, because its 
relaxation time is very long. 

There exist some theoretical studies of coupled oscillators based on the assumption that 
the effect of the frequency disorder is equivalent to that of thermal noise |9|, [10| . According 
to the present results, we can conclude that the effect of quenched frequency disorder is far 
different from that of thermal noise [|ll], [1^ . 
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